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A general algorithm is presented for estimating the nonlinear instability threshold, a c , for subcrit- 
ical transitions in systems where the linearized dynamics is significantly non-normal (i.e. subcritical 
' bifurcations of Takens-Bogdanov type). The JV-dimensional degenerate node is presented as an 

D , example. The predictions are then compared to numerical studies with excellent agreement. 
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PACS numbers: 47.20.Ky, 47.20.Ft, 47.27.Cn 

Consider a nonlinear dynamical system whose dynamics in the vicinity of a stable equilibrium is non-normal at 
linear order. (A matrix or linear operator, C, is non-normal if C^C ^ CC^ with £ T the adjoint.) In the subcritical 
case, nonlincarity makes the equilibrium unstable to finite-amplitude perturbations. The goal is to estimate the size 
of the smallest impulse needed to drive the system unstable, denoted a c . Non- normality implies that the eigenvectors 
of C are not orthogonal. In the extreme case of degeneracy two, or more, of the eigenvectors can become parallel 
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leading to non-diagonalizability of the operator (see, for example, [jl| for a discussion of the finite dimensional case) . 
A prototypical example in two dimensions is that studied by Takens 0] and Bogdanov j^](here the i\ — alx\/dt etc.) 
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with fi and /2 nonlinear functions of their arguments. In a slight abuse of terminology, in the following discussion we 
will use the terms 'Takens-Bogdanov bifurcations' and 'non-normal transitions' interchangeably. Hence, a 'Takens- 
Bogdanov bifurcation' or a 'non-normal transition' here shall mean any system where the linearized dynamics is 
non-normal and the non-normality produces important effects. 

The notion that non-normality might be physically important was first suggested by Orr [|J in an attempt to 
£S1 . explain the failure of standard linear stability analysis to predict the observed critical Reynold's number for the 
laminar/ turbulent transition in some shear flows. This conjecture has been revived recently (see for pro and @ 
for con). Such degenerate transitions occur in other physical models as well. For example they arise in aeroelastic 
models |ll| , and stall models for turbines [^2) . For a discussion of other physical applications the interested reader is 
referred to Chapter 7 of the most recent edition of jl3| . The related degenerate Hopf is discussed in [Q . 

Much of the work cited above deals with bifurcations of low-dimensional models derived from more primitive 
equations (for example Navicr-Stokes) by Galerkin projection. When such projections exhibit degeneracies one must 
be careful about the physical interpretation. Unless one can demonstrate that the behavior of interest is robust under 
perturbation then it is most likely that the degeneracy is a mathematical pathology with little physical importance. 

Such perturbations come in many varieties, with two of most important being: a) perturbations of the equations 
themselves (meaning, for example, the degeneracy is not exact), and b) noise driving, both additive and multiplicative. 
Such noise might represent, for example, coupling to the environment or to degrees of freedom which have been 
projected out by the Galerkin procedure. In this letter we consider perturbations of type a). The additive noise 
response of Takens-Bogdanov systems will be discussed elsewhere [H| . 

A key question from a physical point of view is: What robust observable characteristics, if any, distinguish 'non- 
normal' transitions from 'normal' ones? 

Baggett and Trefethen |6| have shown numerically that a range of low-dimensional non-normal systems exhibit 
anomalous scaling behavior at subcritical transitions: the rate at which a c J, as the threshold for linear instability is 
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approached differs markedly from that of normal systems. This is summarized by the scaling law o~ c ~ e 7 with e the 
linear stability parameter. (Here e is used instead of the inverse Reynolds' number, 1/R. The threshold is e = with 
e > stable.) For normal transitions 7 is generically unity, while for non- normal transitions 7 can be greater than 
unity. Baggett has also shown that it is possible to derive the anomalous scaling exponent in a simple 2-dimcnsional 
case [M. 

The purpose of this letter is to provide a general algorithm for computing the scaling exponents and to illustrate the 
geometric origin of the anomalous scaling behavior. The basic logic parallels that of threshold estimation for normal 
subcritical systems: 

First, at the linear instability threshold (e = 0), a normal form analysis fll3| identifies nonlinear terms which cannot 
be removed - i.e. transformed to higher nonlinear order - by coordinate transformations. These are called resonant 
nonlinear terms. Here, the linear term is assumed to be exactly degenerate and non-diagonalizable. 

Second, backing away from the linear threshold (e > 0), an asymptotic analysis is performed on the resonant 
nonlinear terms to determine which of them dominate. Here, the degeneracy could be weakly broken as well. The 
physical intuition is that the most important nonlinear terms are both resonant (in the normal form) and dominant 
(in the asymptotic limit of e | 0). As will be shown, for the 2-dimcnsional degenerate node this combined normal 
form/asymptotic balance method identifies the same dominant nonlinear term as reported by Takens Takens, 
however, used the technique of topological blowup to identify the dominant nonlinear term, requiring three blowups in 



sequence before the dominant term was revealed |13 



While the method reported here is quite general, it is illustrated on an iV-dimensional degenerate node because this 
is relatively simple. The extension to general N is non-trivial and is necessary for physical applications where more 
than 2 degrees of freedom may be near threshold. To our knowledge, the topological blowup analysis of Takens has 
not been extended to higher dimensions. 

Before starting the analysis, we state our main conclusion: the anomalous scaling behavior of non-normal transitions 
is determined purely by an appropriate balance between the dominant linear and nonlinear effects, as identified by 
the normal form/asymptotic balance analysis. The new element in the non-normal balance is a geometric relation, 
summarized by the triangle relation, Figure (1). We turn now to the presentation of the algorithm. We consider the 
discrete-time map as well as the continuous-time flow to illustrate the simplicity of the result. The flow is gotten by 
an appropriate limit of the map. 




FIG. 1. The asymptotic scaling relation between x\ a and ct c is summarized by the triangle relation of this figure. The 
triangle OS A is formed by the node (O), the saddle (S) and the point of closest approach of the basin boundary to the node 
(which lies in the immediate neighborhood of (A)). The shape of this triangle is determined by the linear amplification factor 
(|): Xi s /a c : \xi max /o~o. 

Consider the discrete-time dynamical system: 

Xj(m + 1) =Fj(x(m);p) to = 0,1,2... (2) 

where i.FeS^f is a smooth nonlinear function of x and p is a control parameter. Suppose F has a fixed point 
%*(p) which we take to be the origin. Expanding to first order: 

N 

x j (m+l) = J2 A jkX k (m); j = l,2...N, (3) 

k=l 
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with A = V.FI 



We take A to be of the form: 
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(4) 



All off-diagonal terms (except a) are zero. The eigenvalue A is real and < A < 1. The coupling constant, a, is 
assumed to be real and positive. More general A will be treated in a longer paper. 

The map (||) can be considered an Euler integrator for the flow x = (A — I)x if we take A — 1 ~ St and a ~ St, with 
St the step-size. Hence, the discrete-time node will have the same scalings as the related flow. In fact, although the 
scaling exponent ( |l4[ ) is developed using maps, the numerical tests (Figure (2)) were all performed using flows. The 
results are in complete agreement. 




FIG. 2. The scaling exponents 7(iV, n) are tested for various combinations of N (the dimension of the dynamics) and n (the 
leading order nonlinearity) . The symbols are the numerical results which were generated using a Bulirsch-Stoer method (ItJ 
and the lines are the predicted scalings. The results are labeled as {N, n, 7). 

Consider the linear impulse response. The analysis for the flow is straightforward as one is dealing with an N- 
dimensional system of equations with constant coefficients. Hence x(t) = exp(At)xQ. The degeneracy of A implies 
that the components of x(t) will not be simple exponents in t, but of the form i n e~ 7t for some 7. This result can 
be found in any undergraduate text on differential equations. The calculation for the discrete-time map is more 
challenging: 

The system is given a random initial kick at m = with correlation matrix: < xoXq >= a^I and I the N x N 
identity matrix. At the mth time step x(m) = A m xo. If A is banded, then so is A m , hence one need only calculate 
the last column, [A m ]j]y: 

[^b„ = (/V) A- (?)-', ,5, 

(m \ 
jy j J the binary coefficients. The index j ranges from N — m to N for m < N (with all entries above the 

j = N — m entry still zero), and from 1 to N when m > N. 

The norm ||a;(m)|| 2 = x(m) ■x(m) — xq ■ [A m ] T A m xo- Taking the ensemble average over the initial conditions gives: 
< ||x(m)|| 2 >= (JqTt (\A T ] m A mS ). For m » N — j the trace of L4 T ]"M m is dominated by the contribution from 
L4 m ]iAr, hence we take that as our estimate: 

< ||x(m)|| 2 > 1 / 2 ^ \ Xl (m)\ ~ a X m ( ^ x ) (^J)^ (6) 

The above result holds for arbitrary A and a as long as a > 0(e). Now assume A = 1 — e and consider e j. 0. Taking 
m ~ t/e with t ~ 0(1) and using Stirling's formula ml ~ TO m e _m (27rm) 1 / 2 , one has 
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As a function of t, this reaches a maximum when t = N — 1, therefore 



cr \6 



(8) 



where we have suppressed the ./V-dependent prefactor. This result is the amplification factor of the linear transients 
(note that it is trivially valid if N = 1). 

Now consider the effects of nonlinearity, and treat first the N = 2 case for simplicity. We assume only that 
the nonlinear terms are smooth. To study generic behavior one casts the given problem into its simplest form by 
performing smooth and invertible coordinate transformations to eliminate as many nonlinear terms as possible. The 
normal form analysis identifies those terms which cannot be eliminated. This analysis is done for e = 0, then e > 
reintroduced for the asymptotic balance estimates to follow. For a detailed discussion of the Takens-Bogdanov normal 
form analysis, the interested reader is refered to Chapter 7 of [ Q . 

To quadratic order, the normal form for the N — 2 degenerate node is 

x[ = (1 - e)xi + ax 2 

x' 2 = (1 - e)x 2 + a 2 xix 2 + b 2 x\ 

with a 2 and b 2 arbitrary coefficients. If these quadratic terms do not appear, then one must go to higher order (all 
other quadratic terms can be pushed to higher order by changing coordinates). At n order one finds 

x[ = (1 - e)xi + ax 2 ^ (1Q . 
x' 2 = (1- e)x 2 +a n x"~ 1 x 2 + b n X^. K ' 

For concreteness assume a n and b n are positive (or 0). This insures the system is subcritical. Consider n = 2: 
solving (^J) for the position of the saddle (find the second root of x' = x): ex± s = ax 2s ; ex 2s = a 2 xi s x 2s + b 2 x\ s . 
There are two simple cases: 

(I) : a 2 = 0,& 2 ~ 0(1); xis = e 2 /ab 2 , x 2s = (e/a)xu- 

(II) : a 2 ~ 0(1), b 2 = 0; x u = e/2a 2 , x 2s = 0. 

Rescaling via x± 8 = e 2 x\ s and x 2s = e 3 x 2s reveals that the b 2 terms dominates. Hence, the scaling for Case I should 
be seen most often for quadratically nonlinear systems, while that of Case II requires special conditions (the smallness 
of b 2 to at least 0(e)). (N.B. As mentioned earlier, this identification of b 2 as the dominant term agrees with the 
topological blowup analysis of Takens p||l5|] . The approach we describe in this letter differs from a blowup analysis in 
that we are considering asymptotic balances in the neighborhood of the bifurcation (e > 0) while the blowup analysis 
is done at the bifurcation (e = 0).) 

At cubic order, it is the 630;^ term driving x 2 which dominates. This behavior holds for general ni/i-order non- 
linearities, which allows us to state: if the first nonlinearities appear at order n, and b n > and ~ 0(1), then the 

position of the saddle is given by ex\ s ~ ax 2s ; ex 2s ~ b n x™ s . This implies Xu ~ ( ^f - ) I x 2s ~ (~f ) ^" s - 

In higher dimensions (N > 2), the normal form analysis reveals that new resonances become possible (with one new 
resonance appearing for each increment N — > TV + 1). Most importantly, the b n x\ term will always resonantly drive 
xn- Asymptotic estimates show that if b n > and ~ 0(1) this term will be dominant. The position of the saddle is: 

exi s - ax 2s ; . . . ex kx ~ ax k +\ ; . • . ex Ns ~ 6„a:™ s , (11) 

which gives 

_e_ /e_\ N ~' 
b n \aJ 



Xls 



l/(n-l) 

(12) 



Note that this estimate is also valid in the normal case where, because the linear term is diagonalizable, N is effectively 
1. 

We now turn to the estimate of the subcritical threshold: as e j 0, how far is the basin boundary from the degenerate 
node? Eq.dH) gives the distance to the basin boundary along the stable manifold of the node. If this were a normal 
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saddle-node bifurcation, x\ s would typically give a good estimate of the distance of closest approach of the basin 
boundary. However, the non-normal linear behavior forces the basin boundary to form an acute angle with the stable 
manifold of the node, hence it will lie very close to the node in directions transverse to the stable manifold. This is 
summarized by the triangle relation of Figure (0). The linear response determines the shape of the triangle and relates 
l^ilmax to an initial impulse x m ~ er via Eq. (||), i.e. |xi| mQa; ~ (a/e) w_1 cr . If |xi| ma:c ~ |a;|i s , then this initial 
perturbation will have crossed the basin boundary. Therefore, this balance determines the threshold for instabilities 
due to finite perturbations and gives as a threshold estimate: 

a c ~ (e/a^-^u = (13) 

with 

, n(N -1) + 1 , s 

7 (N,n) = ^ 14 

n — 1 

Eq.(|l4|) is our primary result. We note that this is also valid for normal (N = 1) case. 

Table [| summarizes the scaling exponents, 7, for several N and n. These were tested numerically with flows. The 
flows have linear dynamics x = Ax with A of the form (^J) and A replaced by — e. The models all have the dominant 
nonlinear term driving the xn component. All nonlinear coefficients were set to unity. The initial conditions were 
set to be (0, • • • ,x° N ). The scaling exponents were computed by plotting the critical x° N beyond which trajectories 
escape to ||a;|| >> \xi s \. The results are summarized in Figure (2). As can be seen, the observed scalings agree 
completely with the prediction (|l4|). 

In summary, we have shown that it is possible to systematically evaluate the importance of various nonlinear effects 
on non-normal transitional behavior by an extension of the techniques used for normal systems. This leads to an 
algorithm capable of predicting the nonlinear threshold for subcritical transitions. The algorithm was illustrated by 
application to an JV-dimensional degenerate node, where the anomalous scaling behavior was shown to be due to the 
fact that the non-normality of the linear term forces a geometrical relationship between the length scales along and 
across the stable manifold of the node, an effect which is absent in normal systems. The resulting scaling exponent (|l4| ) 
shows that non-normality and nonlinearity act together to increase the sensitivity to subcritical transitions, and the 
threshold depends exponentially on the number of degrees of freedom taking part, N. Figure (g) suggests that by 
measuring such scaling behavior near threshold it might be possible to choose between various models (or at least 
eliminate a large class of them), though if N is large such scaling regimes will be extremely narrow in e. 

We thank J. S. Baggett for useful comments. This work was supported by the AFOSR and the DOE. 
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TABLE I. Tabular summary of 7 computed from ([h]). The highlighted entries correspond to models in ||, with (their 
notation) *=TTRD', **=TTRD". As mentioned in the text, although both of their models are nominally quadratic, the 
normal forms are quite different and show that TTRD" can in fact be transformed to be cubically nonlinear. Their reported 
threshold scalings for TTRD' and TTRD" are 3 and 2, respectively. These are the only two models we compare with Q because 
the rest of their models either have a square root singularity at the origin - hence the normal form analysis does not apply - 
or the models are not uniform in their coupling, implying the normal form used here would not be the correct one. 
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